BARRETO et al. 2022: Habitat loss and biodiversity gain (BIOCON-S-22-00746)

Ver com a Me:

  1. conferência dos recálculos, alpha FS parou de responder
  2. Olhar o código e resíduos, alguns já não parecem mais ok
  3. marcar na planilha os bichos que mudaram
  4. cálculo beta RC dentro do wrangling ou separado? onde jogo as funções (em script separado e usar source)
  5. sugestão de como plotar quando tem nulo (linha tracejada mais clara e sem SE?)
  6. Plot também vai pro Git? script separado?

Model formulas

For each of our response variables we run three candidate models for each scale:

Model 0 <- ~1, absence of effect

Model 1 <- percentage of forest cover at 3 and 5km radius

Model1.2 <- non-linear forest cover at 3 and 5km radius (quadratic model)

Loading data and gathering into one dataset
Landscape gamma_FS gamma_NFS alpha_FS alpha_NFS betad_FS betad_NFS ab.land_FS ab.land_NFS brc.fsNA brc.nfsNA brc.abund.fsNA brc.abund.nfsNA fc_3km fc_5km
148 11 9 3.75 2.25 7.25 6.75 66 40 -0.2928571 0.1826667 -0.0014658 0.4252252 38.52523 43.86337
215 19 7 5.62 2.00 13.38 5.00 174 51 0.1057143 -0.3123810 0.2759545 0.3903427 30.15485 25.23160
263 16 13 6.00 4.50 10.00 8.50 393 214 -0.4919048 -0.0323810 0.6333476 0.3923447 10.10641 13.41180
266 11 8 3.25 1.38 7.75 6.62 70 12 -0.2257143 0.1338095 0.1087754 0.1472425 30.87525 31.16090
282 17 11 6.50 4.12 10.50 6.88 219 212 -0.3275000 -0.1152381 0.2126412 0.7177654 15.26270 16.35308
291 9 3 3.71 1.71 5.29 1.29 150 49 -0.3866667 -0.5904762 0.1027218 -0.4251871 45.74294 45.13213
317 12 4 5.00 1.38 7.00 2.62 152 56 -0.3853571 -0.5321429 0.1819677 -0.4317889 48.75805 47.85138
329 21 14 8.62 5.75 12.38 8.25 369 262 -0.0871429 0.2217857 0.9169527 0.4891320 14.73791 13.62588
333 14 15 4.50 4.62 9.50 10.38 158 219 -0.1275000 0.1932143 0.6788932 0.6190476 13.71403 14.50960
335 15 8 5.75 2.38 9.25 5.62 140 113 -0.2564286 -0.3366667 0.3479193 -0.9437533 23.87913 19.21670
359 16 17 6.00 5.50 10.00 11.50 149 133 -0.2253571 0.3128571 0.7704490 0.8121693 18.10208 17.62431
399 20 14 7.88 5.12 12.12 8.88 432 169 -0.1223810 0.2703571 0.9618190 0.6938009 29.73616 24.97851

RESULTS

We found 28 species to be forest specialists, those that come from the Atlantic Forest domain and exclusively collected in forested habitats; while 23 were non-forest specialists, species from a broader domain and/or less restrictive to forest habitats (hereafter FS and NFS, respectively).

FOREST SPECIALISTS

Gamma diversity

Modelling

g0 <- glm(gamma_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(g.fs.aic <- AICctab(g0, g5k, g5kq, g3k, g3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## g5k  -27.9   64.8   4.7     0.0 3  0.731 
## g3k  -29.6   68.2   3.0     3.4 3  0.136 
## g5kq -27.8   69.4   4.8     4.6 4  0.074 
## g0   -32.6   70.6   0.0     5.8 2  0.041 
## g3kq -29.2   72.2   3.4     7.3 4  0.019

Summary top model

summary(g5k)
## 
## Call:
## glm(formula = gamma_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.23181  -0.12743  -0.03948   0.15159   0.27856  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 20.68384    1.96964   10.50 1.01e-06 ***
## fc_5km      -0.21451    0.05782   -3.71  0.00404 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.0352002)
## 
##     Null deviance: 0.74994  on 11  degrees of freedom
## Residual deviance: 0.34328  on 10  degrees of freedom
## AIC: 61.839
## 
## Number of Fisher Scoring iterations: 3

Inspect residuals:

sim <- simulateResiduals(g5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1487, p-value = 0.632
## alternative hypothesis: two.sided
hnp::hnp(g5k)
## Gamma model

boot::glm.diag.plots(g5k) ## for glm poisson or neg binom

Prediction:

Alpha diversity

Modelling

a0 <- glm(alpha_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(a.fs.aic <- AICctab(a0, a5k, a5kq,
                     a3k, a3kq, 
                     base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## a5k  -19.3   47.6   2.8     0.0 3  0.523 
## a3k  -20.3   49.5   1.9     1.9 3  0.199 
## a0   -22.1   49.6   0.0     2.0 2  0.195 
## a5kq -19.0   51.8   3.1     4.2 4  0.064 
## a3kq -20.3   54.2   1.9     6.6 4  0.019

Summary top model

summary(a5k)
## 
## Call:
## glm(formula = alpha_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.43224  -0.11774  -0.03674   0.08987   0.35754  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.42155    0.94502   7.853 1.38e-05 ***
## fc_5km      -0.07200    0.02815  -2.558   0.0285 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06034959)
## 
##     Null deviance: 0.97549  on 11  degrees of freedom
## Residual deviance: 0.61269  on 10  degrees of freedom
## AIC: 44.608
## 
## Number of Fisher Scoring iterations: 4

Inspect residuals:

sim <- simulateResiduals(a5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1043, p-value = 0.656
## alternative hypothesis: two.sided
hnp::hnp(a5k)
## Gamma model

boot::glm.diag.plots(a5k)

Prediction:

Beta diversity additive

Modelling

### LM Beta for regular model selection
b0 <- glm(betad_FS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(b.fs.aic <- AICctab(b0, b5k, b5kq, 
                     b3k, b3kq, 
                     base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## b5k  -22.1   53.3   5.1     0.0 3  0.690 
## b5kq -21.4   56.6   5.8     3.3 4  0.133 
## b3k  -23.9   56.8   3.3     3.6 3  0.116 
## b3kq -22.7   59.2   4.5     5.9 4  0.036 
## b0   -27.2   59.8   0.0     6.6 2  0.026

Summary top model

summary(b5k)
## 
## Call:
## glm(formula = betad_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.24293  -0.12664  -0.05680   0.08772   0.34333  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.29031    1.24634  10.663  8.8e-07 ***
## fc_5km      -0.14349    0.03623  -3.961  0.00268 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.03508083)
## 
##     Null deviance: 0.76409  on 11  degrees of freedom
## Residual deviance: 0.32769  on 10  degrees of freedom
## AIC: 50.257
## 
## Number of Fisher Scoring iterations: 4

Inspect residuals:

sim <- simulateResiduals(b5k)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1585, p-value = 0.608
## alternative hypothesis: two.sided
hnp::hnp(b5k) ## for glm Gamma or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(b5k) ## for glm Gamma or neg binom

Prediction

bRaup-Crick (presence/absence-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.fsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.fsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq, 
                      rc3k, rc3kq,
                      base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc0    5.4   -5.4  0.0     0.0  2  0.475 
## rc3kq  8.7   -3.7  3.3     1.7  4  0.203 
## rc5k   5.9   -2.7  0.5     2.7  3  0.124 
## rc5kq  8.1   -2.6  2.8     2.9  4  0.113 
## rc3k   5.5   -2.0  0.1     3.4  3  0.085

Summary of 2nd best model, as the top is of absence of effect (~1)

summary(rc3kq)
## 
## Call:
## lm(formula = brc.fsNA ~ fc_3km + I(fc_3km^2), data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.13500 -0.09118 -0.03185  0.07770  0.23076 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.6988850  0.2220491  -3.147   0.0118 *
## fc_3km       0.0413030  0.0173151   2.385   0.0409 *
## I(fc_3km^2) -0.0007386  0.0002921  -2.529   0.0323 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.135 on 9 degrees of freedom
## Multiple R-squared:  0.4268, Adjusted R-squared:  0.2995 
## F-statistic: 3.351 on 2 and 9 DF,  p-value: 0.08171

Inspect residuals:

sim <- simulateResiduals(rc3kq)
plot(sim)
## qu = 0.25, log(sigma) = -3.759907 : outer Newton did not converge fully.

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.82544, p-value = 0.816
## alternative hypothesis: two.sided
hnp::hnp(rc3kq)
## Gaussian model (lm object)

Prediction:

Total abundance

Modelling

### LM "gamma"  abundanec for regular model selection
gab0 <- glm.nb(ab.land_FS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_FS ~ fc_5km +  I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_FS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_FS ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(gab.fs.aic <- AICctab(gab0, gab5k, gab5kq, 
                       gab3k, gab3kq,
                       base=T,weights=T, logLik=T))
##        logLik AICc  dLogLik dAICc df weight
## gab5k  -70.5  150.1   2.1     0.0 3  0.399 
## gab0   -72.6  150.5   0.0     0.5 2  0.318 
## gab3k  -71.2  151.3   1.4     1.2 3  0.217 
## gab5kq -70.4  154.5   2.2     4.5 4  0.043 
## gab3kq -71.0  155.8   1.6     5.7 4  0.023

Summary 2nd best model, as top model was the one of absence of effect (~1):

summary(gab5k)
## 
## Call:
## glm.nb(formula = ab.land_FS ~ fc_5km, data = data_div, init.theta = 4.537178159, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6496  -0.9873  -0.2873   0.6438   1.8333  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.93023    0.31931  18.572   <2e-16 ***
## fc_5km      -0.02485    0.01110  -2.239   0.0251 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(4.5372) family taken to be 1)
## 
##     Null deviance: 17.268  on 11  degrees of freedom
## Residual deviance: 12.411  on 10  degrees of freedom
## AIC: 147.08
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  4.54 
##           Std. Err.:  1.83 
## 
##  2 x log-likelihood:  -141.082

Inspect residuals:

sim <- simulateResiduals(gab5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1614, p-value = 0.56
## alternative hypothesis: two.sided
hnp::hnp(gab5k)
## Negative binomial model (using MASS package)

boot::glm.diag.plots(gab5k)

# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab3k.r2 <- r2(gab3k)

## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance FS", 
                            "FC at 5km", 
                            round(gab.fs.aic$AICc[1],2),
                            round(gab.fs.aic$dAICc[1],1), 
                            gab.fs.aic$df[1],
                            round(gab.fs.aic$weight[1],2),
                            round(gab5k.r2$R2_Nagelkerke,2)),
                 c(" ", "~ 1 (NULL)", 
                   round(gab.fs.aic$AICc[2],2), 
                   round(gab.fs.aic$dAICc[2],1),
                   gab.fs.aic$df[2], 
                   round(gab.fs.aic$weight[2],2),
                    "-"),
                 c(" ", "FC at 3km", 
                            round(gab.fs.aic$AICc[3],2),
                            round(gab.fs.aic$dAICc[3],1), 
                            gab.fs.aic$df[3],
                            round(gab.fs.aic$weight[3],2),
                            round(gab3k.r2$R2_Nagelkerke,2)))

Prediction:

bRaup-Crick (abundance-based)

Modelling:

### LM gamma for regular model selection
rc0 <- lm(brc.abund.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.fsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.fsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq, 
                      rc3k, rc3kq,
                      base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc5k   0.2    8.6  3.8     0.0  3  0.593 
## rc3k  -0.7   10.4  2.9     1.8  3  0.237 
## rc0   -3.6   12.5  0.0     3.9  2  0.085 
## rc5kq  0.3   13.1  3.9     4.5  4  0.062 
## rc3kq -0.7   15.2  2.9     6.5  4  0.022

Summary top model:

summary(rc5k)
## 
## Call:
## lm(formula = brc.abund.fsNA ~ fc_5km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.39458 -0.18080 -0.00713  0.15192  0.50954 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.900965   0.175354   5.138 0.000439 ***
## fc_5km      -0.017963   0.006073  -2.958 0.014339 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2607 on 10 degrees of freedom
## Multiple R-squared:  0.4666, Adjusted R-squared:  0.4133 
## F-statistic: 8.749 on 1 and 10 DF,  p-value: 0.01434

Inspect residuals:

sim <- simulateResiduals(rc5k)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc5k)
## Gaussian model (lm object)

Prediction

NON-FOREST SPECIALISTS

Gamma diversity

Modelling:

### LM alpha for regular model selection
g0 <- glm(gamma_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(g.nfs.aic <- AICctab(g0, 
                      g5k, g5kq, 
                      g3k, g3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## g3k  -27.3   63.5   7.6     0.0 3  0.7911
## g5k  -29.2   67.5   5.7     4.0 3  0.1095
## g3kq -27.1   68.0   7.8     4.4 4  0.0864
## g5kq -29.2   72.2   5.7     8.6 4  0.0105
## g0   -34.9   75.2   0.0    11.6 2  0.0024

Summary top model:

summary(g3k)
## 
## Call:
## glm(formula = gamma_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.38392  -0.22774  -0.05573   0.15121   0.43058  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.37185    2.13069   8.622 6.07e-06 ***
## fc_3km      -0.30280    0.05227  -5.793 0.000175 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.08050867)
## 
##     Null deviance: 2.71637  on 11  degrees of freedom
## Residual deviance: 0.78116  on 10  degrees of freedom
## AIC: 60.54
## 
## Number of Fisher Scoring iterations: 4

Inspect residuals:

sim <- simulateResiduals(g3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90839, p-value = 0.96
## alternative hypothesis: two.sided
hnp::hnp(g3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(g3k) ## for glm poisson or neg binom

Prediction:

Alpha diversity

Modelling:

### LM alpha for regular model selection
a0 <- glm(alpha_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(a.nfs.aic <- AICctab(a0, a5k, a5kq, 
                      a3k, a3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## a3k  -16.7   42.3   5.7     0.0 3  0.490 
## a5k  -16.9   42.9   5.4     0.6 3  0.371 
## a5kq -16.2   46.2   6.1     3.9 4  0.071 
## a3kq -16.5   46.6   5.9     4.3 4  0.057 
## a0   -22.3   50.0   0.0     7.7 2  0.011

Summary top model:

summary(a3k)
## 
## Call:
## glm(formula = alpha_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.68238  -0.15994   0.00108   0.11140   0.55352  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.82215    0.83950   6.935 4.02e-05 ***
## fc_3km      -0.09207    0.02113  -4.357  0.00143 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1163241)
## 
##     Null deviance: 3.1391  on 11  degrees of freedom
## Residual deviance: 1.2537  on 10  degrees of freedom
## AIC: 39.336
## 
## Number of Fisher Scoring iterations: 6

Inspect residuals:

sim <- simulateResiduals(a3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.93258, p-value = 0.928
## alternative hypothesis: two.sided
hnp::hnp(a3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(a3k) ## for glm poisson or neg binom

Prediction:

Beta diversity additive

Modelling

### LM Beta for regular model selection
b0 <- glm(betad_NFS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(b.nfs.aic <- AICctab(b0, b5k, b5kq, 
                      b3k, b3kq, 
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## b3k  -25.0   58.9   5.9     0.0 3  0.753 
## b5k  -26.8   62.7   4.0     3.7 3  0.116 
## b3kq -24.6   62.9   6.3     3.9 4  0.105 
## b5kq -26.6   67.0   4.2     8.1 4  0.013 
## b0   -30.9   67.1   0.0     8.2 2  0.013

Summary top model:

summary(b3k)
## 
## Call:
## glm(formula = betad_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.71147  -0.22545  -0.03265   0.18829   0.45483  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.59989    1.68741   7.467 2.14e-05 ***
## fc_3km      -0.21226    0.04078  -5.205 0.000399 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1116167)
## 
##     Null deviance: 3.2119  on 11  degrees of freedom
## Residual deviance: 1.2312  on 10  degrees of freedom
## AIC: 55.933
## 
## Number of Fisher Scoring iterations: 4

Inspect residuals:

sim <- simulateResiduals(b3k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.65877, p-value = 0.68
## alternative hypothesis: two.sided
hnp::hnp(b3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(b3k) ## for glm poisson or neg binom

Prediction:

bRaup-Crick (presence/absence-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.nfsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.nfsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq, 
                       rc3k, rc3kq,
                       base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc3k  -0.6   10.2  2.3     0.0  3  0.390 
## rc0   -2.9   11.1  0.0     0.9  2  0.250 
## rc5k  -1.2   11.4  1.7     1.2  3  0.219 
## rc3kq  0.5   12.7  3.4     2.5  4  0.114 
## rc5kq -0.9   15.5  2.0     5.3  4  0.027

Summary of 2nd best model, 1st wqs of absence of effect (~1):

summary(rc3k)
## 
## Call:
## lm(formula = brc.nfsNA ~ fc_3km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32488 -0.21626 -0.05457  0.24366  0.39972 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.322872   0.191341   1.687   0.1224  
## fc_3km      -0.014015   0.006519  -2.150   0.0571 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2787 on 10 degrees of freedom
## Multiple R-squared:  0.3161, Adjusted R-squared:  0.2477 
## F-statistic: 4.622 on 1 and 10 DF,  p-value: 0.05707

Inspect residuals:

sim <- simulateResiduals(rc3k)
plot(sim)
## qu = 0.75, log(sigma) = -3.658389 : outer Newton did not converge fully.

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)
## Gaussian model (lm object)

Prediction:

Total abundance

Modelling:

### LM "gamma"  abundanec for regular model selection
gab0 <- glm.nb(ab.land_NFS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_NFS ~ fc_5km +  I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_NFS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_NFS ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(gab.nfs.aic <- AICctab(gab0, gab5k, gab5kq, 
                        gab3k, gab3kq,
                        base=T,weights=T, logLik=T))
##        logLik AICc  dLogLik dAICc df weight
## gab5k  -63.9  136.9   5.3     0.0 3  0.404 
## gab3k  -64.3  137.6   4.9     0.7 3  0.284 
## gab5kq -62.1  137.8   7.2     0.9 4  0.252 
## gab3kq -63.7  141.2   5.5     4.3 4  0.047 
## gab0   -69.2  143.8   0.0     6.9 2  0.013

Summary top model:

summary(gab5k)
## 
## Call:
## glm.nb(formula = ab.land_NFS ~ fc_5km, data = data_div, init.theta = 4.164123767, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8516  -0.4183   0.2175   0.4846   0.8634  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.93704    0.33661  17.638  < 2e-16 ***
## fc_5km      -0.04838    0.01182  -4.094 4.24e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(4.1641) family taken to be 1)
## 
##     Null deviance: 29.266  on 11  degrees of freedom
## Residual deviance: 12.765  on 10  degrees of freedom
## AIC: 133.88
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  4.16 
##           Std. Err.:  1.77 
## 
##  2 x log-likelihood:  -127.88

Inspect residuals:

sim <- simulateResiduals(gab5k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.45825, p-value = 0.288
## alternative hypothesis: two.sided
hnp::hnp(gab5k)
## Negative binomial model (using MASS package)

boot::glm.diag.plots(gab5k)

# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab5kq.r2 <- r2(gab5kq)
gab3k.r2 <- r2(gab3k)

## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance NFS", 
                            "FC at 5km", 
                            round(gab.nfs.aic$AICc[1],2),
                            round(gab.nfs.aic$dAICc[1],1), 
                            gab.nfs.aic$df[1],
                            round(gab.nfs.aic$weight[1],2),
                            round(gab5k.r2$R2_Nagelkerke,2)),
                 c(" ", "FC at 3km",
                   round(gab.nfs.aic$AICc[2],2), 
                   round(gab.nfs.aic$dAICc[2],1),
                   gab.nfs.aic$df[2], 
                   round(gab.nfs.aic$weight[2],2),
                   round(gab3k.r2$R2_Nagelkerke,2)),
                 c(" ", "Non-linear FC at 5km",
                   round(gab.nfs.aic$AICc[3],2), 
                   round(gab.nfs.aic$dAICc[3],1),
                   gab.nfs.aic$df[3], 
                   round(gab.nfs.aic$weight[3],2),
                   round(gab5kq.r2$R2_Nagelkerke,2)
                 ))

Prediction:

bRaup-Crick (abundance-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.abund.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.nfsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.nfsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq, 
                       rc3k, rc3kq,
                       base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc3k  -7.4   23.8  2.0     0.0  3  0.381 
## rc0   -9.4   24.1  0.0     0.2  2  0.337 
## rc5k  -8.0   24.9  1.4     1.1  3  0.219 
## rc3kq -7.3   28.3  2.1     4.5  4  0.041 
## rc5kq -7.9   29.5  1.5     5.7  4  0.022

Summary 2nd best model, 1st is of absence of effect (~1):

summary(rc3k)
## 
## Call:
## lm(formula = brc.abund.nfsNA ~ fc_3km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.24646 -0.18490  0.04466  0.26676  0.52334 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.84188    0.33775   2.493   0.0318 *
## fc_3km      -0.02258    0.01151  -1.962   0.0781 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4919 on 10 degrees of freedom
## Multiple R-squared:  0.278,  Adjusted R-squared:  0.2058 
## F-statistic: 3.851 on 1 and 10 DF,  p-value: 0.07814

Inspect residuals:

sim <- simulateResiduals(rc3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)
## Gaussian model (lm object)

Prediction:

AICc-model selection table

Table S2. AICc-based selected models table containing the models of gamma (γ), alpha (ɑ), additive beta (β) diversities followed by abundance at landscape total, abundance at mean site and β-diversity calculated from the Raup-Crick null model approach, for both presence/absence (βRC) and abundance-based (βRC-abundance) that scored below ΔAICc=2. All tested as a function of Atlantic Forest loss for each group of species, first models separated into forest specialists species (FS; Table 1a), and followed by non-forest species (NFS; Table 1b). We show selected models ordered from the lowest to the highest AICc values ΔAICc≤ 2. We present information regarding the diversity metric tested (Diversity), the model terms for forest loss predictors (Model), the Akaike information criterion corrected for small samples (AICc), the AICc difference from the first-ranked model (ΔAICc), degrees of freedom (df) Akaike weight (ωi) and model variance explained (r²).
Diversity Model AICc dAICc df weight r2
Gamma FS FC at 5km 64.84 0 3 0.73 0.55
Alpha FS FC at 5km 47.61 0 3 0.52 0.38
FC at 3km 49.54 1.9 3 0.2 0.27
~ 1 (NULL) 49.58 2 2 0.19
Beta FS FC at 5km 53.26 0 3 0.69 0.58
RCbeta (P/A-based) FS ~ 1 (NULL) -5.44 0 2 0.47
Non-linear FC at 3km -3.74 1.7 4 0.2 0.3
FC at 5km -2.75 2.7 3 0.12 -0.01
Total abundance FS FC at 5km 150.08 0 3 0.4 0.44
~ 1 (NULL) 150.54 0.5 2 0.32
FC at 3km 151.31 1.2 3 0.22 0.32
RCbeta (abund-based) FS FC at 5km 8.61 0 3 0.59 0.41
FC at 3km 10.44 1.8 3 0.24 0.32
Gamma NFS FC at 3km 63.54 0 3 0.79 0.74
Alpha NFS FC at 3km 42.34 0 3 0.49 0.63
FC at 5km 42.89 0.6 3 0.37 0.61
Beta NFS FC at 3km 58.93 0 3 0.75 0.65
RCbeta (P/A-based) NFS FC at 3km 10.2 0 3 0.39 0.17
~ 1 (NULL) 11.09 0.9 2 0.25
FC at 5km 11.35 1.2 3 0.22 0.17
Total abundance NFS FC at 5km 136.88 0 3 0.4 0.82
FC at 3km 137.58 0.7 3 0.28 0.79
Non-linear FC at 5km 137.82 0.9 4 0.25 0.93
RCbeta (abund-based) NFS FC at 3km 23.84 0 3 0.38 0.21
~ 1 (NULL) 24.08 0.2 2 0.34
FC at 5km 24.94 1.1 3 0.22 0.13

FIGURES

Gamma hab

### Gamma ####
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))

## prediction g2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(g5k, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
## prediction g1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds2 <- predict(g3k, type="response", newdata=new2, se.fit=T)

new2$pred <- preds2$fit
new2$ic.up <- preds2$fit + preds2$se.fit*1.96 
new2$ic.low <- preds2$fit - preds2$se.fit*1.96 

## Plot gamma richness ~ fc %
(g.hab <- ggplot(data= new, aes(x= fc_5km, y = pred)) + 
    geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
    geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
    geom_point(data= data_div, aes(x=fc_5km, y=gamma_FS), col= "#0072B2") + 
    geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
    geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
    geom_point(data= data_div, aes(x=fc_3km, y=gamma_NFS), col= "#D55E00") +
    xlab(" ") + ylab("Gamma") + 
    theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
    ylim(0,max(data_div[2:7])+0.5) + 
    # scale_x_reverse() +
    annotate("text", y=20.5, x=38,colour= "#0072B2",
             label=expression(paste(~ R^2 ~ "= 0.66"))) +
    annotate("text", y=19, x=40, colour= "#D55E00",
             label=expression(paste(~ R^2 ~ "= 0.74"))) + labs(tag="a"))

Alpha hab

a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
## prediction a2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(a5k, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
## prediction a1
new2 <-  data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(a3k, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot Alpha richness ~ fc_3km
(a.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=alpha_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=alpha_NFS), col= "#D55E00") +
  xlab("") + ylab("Alpha") + 
  scale_color_manual(breaks=c("FS", "NFS"),
                     values=c('FS'='#0072B2', 'NFS'='#D55E00')) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=14,face="bold"),
        legend.position = 'bottom') + 
  ylim(0,max(data_div[2:7])+0.5) +
    # scale_x_reverse() +
  annotate("text", y=13, x=38,colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0383")))+
  annotate("text", y=11, x=40, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.63"))) +
  labs(tag="b"))

Beta hab

### Beta add  ####
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
## prediction b1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(b5k, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))

## prediction b1
new2 <-  data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(b3k, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot Beta Diversity ~ fc %
(b.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=betad_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x= fc_3km, y= pred), 
            col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=betad_NFS), col= "#D55E00") +
  xlab(" ") + ylab("Beta") + 
  theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
  ylim(0,max(data_div[2:7])+0.5) + 
    # scale_x_reverse() +
  annotate("text", y=18, x=40,colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.58")))+
  annotate("text", y=16, x=38, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.65"))) +
  labs(tag="c"))

## plot here not to need posterior figure editing app
g.hab + a.hab + b.hab

#ggsave("figures/plot.hab.jpeg",units="cm", width=20, height=8, dpi=300)

Total abundance FS & NFS

## FS
gab1 <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
## prediction gab1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(gab1, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
gab1 <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
## prediction gab2
new2 <- new
preds <- predict(gab1, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot gamma richness ~ forest cover 5km
(gab.hab <- ggplot(new, aes(x=fc_5km, y=pred)) + 
  geom_line(linetype= "dashed", col= "#0072B2", aes(alpha= 0.1)) +
  # geom_ribbon(aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1,  fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=ab.land_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_5km, y=pred), col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.2, fill= "#D55E00") +
  geom_point(data= data_div, aes(x=fc_5km, y=ab.land_NFS), col= "#D55E00") +
    scale_x_continuous(limits= c(10,50)) +
  xlab(" ") + ylab("Total abundance") + 
  theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
  # scale_x_reverse() +
    annotate("text", y=300, x=38, alpha= .5, colour= "#0072B2",
             label=expression(paste(~ R^2 ~ "= 0.44"))) +
    annotate("text", y=270, x=40, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.82"))) +
  labs(tag="d"))

### Beta RC FS & NFS

## FS
rc2 <- lm(brc.fsNA ~ fc_3km +  I(fc_3km^2), data=data_div)
## prediction rc2.2 and rc2
new <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=30))
preds <- predict(rc2, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
rc1 <- lm(brc.nfsNA ~ fc_3km, data=data_div)
## prediction rc2
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot RCbeta abundance ~ forest cover at 3
(rc.hab <- ggplot(data= new, aes(x=fc_3km, y=pred)) + 
  geom_line(data= new,aes(x=fc_3km, y=pred), 
            col= "#0072B2", linetype= "dashed", alpha= .2) +
  # geom_ribbon(data= new, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") + 
  geom_point(data=data_div, aes(x=fc_3km, y=brc.fsNA), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D33E00", linetype= "dashed", alpha= .2) +
  # geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D33E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=brc.nfsNA), col= "#D33E00") +
  xlab("Percent forest cover") + ylab(c(expression(bold("β"["RC"])))) + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),
        legend.position = 'none') + ylim(-1,1) + 
    # scale_x_reverse() +
  annotate("text", y=0.9, x=43, alpha= .5, colour= "#D33E00",
           label=expression(paste(~ R^2 ~ "= 0.82"))) +
  annotate("text", y=0.7, x=43, alpha= .5, colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.44"))) +
  # annotate("text", y=0.9, x=43, colour= "#D33E00",
  #          label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.37"))) +
  # annotate("text", y=0.8, x=43, colour= "#0072B2",
  #          label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.33"))) +
  labs(tag="e"))

Beta RC abund FS & NFS

rc2 <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
## prediction rc2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(rc2, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
rc1 <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
## prediction rc1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## RCbeta abund ~ FC%
(rc.ab.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new,aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") + 
  geom_point(data=data_div, aes(x=fc_5km, y=brc.abund.fsNA), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), linetype= "dashed", col= "#D55E00", alpha=.2) +
  # geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=brc.abund.nfsNA), col= "#D55E00") +
  xlab(" ") + ylab(c(expression(bold("β"["RC-abundance"])))) + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
    ylim(-1,1.1) + 
    # scale_x_reverse() +
  annotate("text", y=0.9, x=45, colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.41"))) +
  annotate("text", y=0.75, x=43, alpha= .5, colour= "#D55E00",
           label=expression(paste( ~ R^2 ~ "= 0.21"))) +
  labs(tag="f"))

Final plot

(plot.hab <- (g.hab + a.hab + b.hab) /
    (gab.hab + rc.hab + rc.ab.hab))

#ggsave("figures/figure2.jpeg",units="cm", width=28, height=15, dpi=300)